R Markdown
library(ggplot2)
library(tidyverse)
library(dplyr)
library(meta)
library(PRISMAstatement)
library(skimr)
library(MASS)
library(ggpubr)
library(patchwork)
### Site Level
final_data <- read.csv("Final Data.csv")
final_data_2022 <- final_data %>%
filter(year == "2022")
shapiro.test(final_data_2022$animals)
##
## Shapiro-Wilk normality test
##
## data: final_data_2022$animals
## W = 0.76017, p-value = 7.037e-05
ggqqplot(final_data_2022$animals)
n1 <- glm(animals ~ shrub_density * aridity, family = "gaussian", data = final_data_2022)
n1
##
## Call: glm(formula = animals ~ shrub_density * aridity, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density aridity
## -113.435 -6.222 50.977
## shrub_density:aridity
## 2.798
##
## Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
## Null Deviance: 619900
## Residual Deviance: 517700 AIC: 317.6
anova(n1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 619927
## shrub_density 1 26978 22 592949 0.30732
## aridity 1 71496 21 521453 0.09653 .
## shrub_density:aridity 1 3720 20 517733 0.70462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(n1), "\n")
## Linear Model AIC: 317.6089
quadratic_model_1 <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2022)
quadratic_model_1
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## -4704.343 69.124 -7.773
## aridity I(aridity^2) shrub_density:aridity
## 2395.612 -287.024 6.017
##
## Degrees of Freedom: 23 Total (i.e. Null); 18 Residual
## Null Deviance: 619900
## Residual Deviance: 180300 AIC: 296.3
anova(quadratic_model_1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 619927
## shrub_density 1 26978 22 592949 0.100791
## I(shrub_density^2) 1 2942 21 590007 0.587895
## aridity 1 78131 20 511876 0.005227 **
## I(aridity^2) 1 315016 19 196860 2.052e-08 ***
## shrub_density:aridity 1 16537 18 180323 0.198864
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_1), "\n")
## Quadratic Model AIC: 296.2959
if (AIC(quadratic_model_1) < AIC(n1)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
### Richness
n2 <- glm(richness ~ shrub_density * aridity, family = "gaussian", data = final_data_2022)
n2
##
## Call: glm(formula = richness ~ shrub_density * aridity, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density aridity
## -7.07856 0.25745 3.07545
## shrub_density:aridity
## -0.04175
##
## Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
## Null Deviance: 475
## Residual Deviance: 334.1 AIC: 141.3
anova(n2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 474.96
## shrub_density 1 6.431 22 468.53 0.534951
## aridity 1 133.613 21 334.91 0.004681 **
## shrub_density:aridity 1 0.828 20 334.09 0.823795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(n2), "\n")
## Linear Model AIC: 141.3094
quadratic_model_2 <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2022)
quadratic_model_2
##
## Call: glm(formula = richness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## -128.52243 -0.17675 0.03871
## aridity I(aridity^2) shrub_density:aridity
## 65.21718 -7.60162 -0.05383
##
## Degrees of Freedom: 23 Total (i.e. Null); 18 Residual
## Null Deviance: 475
## Residual Deviance: 26.87 AIC: 84.82
anova(quadratic_model_2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 474.96
## shrub_density 1 6.431 22 468.53 0.03794 *
## I(shrub_density^2) 1 106.430 21 362.10 < 2e-16 ***
## aridity 1 101.019 20 261.08 < 2e-16 ***
## I(aridity^2) 1 232.882 19 28.20 < 2e-16 ***
## shrub_density:aridity 1 1.324 18 26.87 0.34636
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_2), "\n")
## Quadratic Model AIC: 84.82182
if (AIC(quadratic_model_2) < AIC(n2)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
### Evenness
n3 <- glm(Average_Evenness ~ shrub_density * aridity, data = final_data_2022)
n3
##
## Call: glm(formula = Average_Evenness ~ shrub_density * aridity, data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density aridity
## -0.059619 -0.002519 0.035851
## shrub_density:aridity
## 0.001281
##
## Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
## Null Deviance: 0.1455
## Residual Deviance: 0.1074 AIC: -51.7
anova(n3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.14552
## shrub_density 1 0.0073326 22 0.13819 0.24267
## aridity 1 0.0299710 21 0.10822 0.01817 *
## shrub_density:aridity 1 0.0007791 20 0.10744 0.70332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(n3), "\n")
## Linear Model AIC: -51.70484
quadratic_model_3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2022)
quadratic_model_3
##
## Call: glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## -1.8496731 -0.0253566 0.0022251
## aridity I(aridity^2) shrub_density:aridity
## 0.9526101 -0.1121067 0.0004442
##
## Degrees of Freedom: 23 Total (i.e. Null); 18 Residual
## Null Deviance: 0.1455
## Residual Deviance: 0.01326 AIC: -97.91
anova(quadratic_model_3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.145519
## shrub_density 1 0.007333 22 0.138186 0.001608 **
## I(shrub_density^2) 1 0.056383 21 0.081804 < 2.2e-16 ***
## aridity 1 0.018828 20 0.062975 4.308e-07 ***
## I(aridity^2) 1 0.049621 19 0.013354 2.285e-16 ***
## shrub_density:aridity 1 0.000090 18 0.013264 0.726529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_2), "\n")
## Quadratic Model AIC: 84.82182
if (AIC(quadratic_model_2) < AIC(n2)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
even_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_dens_2022
plot <- abund_dens_2022/rich_dens_2022/even_dens_2022
plot
### Site Level
final_data_2023 <- final_data %>%
filter(year == "2023")
x1 <- glm(animals ~ shrub_density * aridity, family = "gaussian", data = final_data_2023)
x1
##
## Call: glm(formula = animals ~ shrub_density * aridity, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density aridity
## 3.96923 -0.08162 0.93546
## shrub_density:aridity
## 0.20460
##
## Degrees of Freedom: 21 Total (i.e. Null); 18 Residual
## Null Deviance: 12650
## Residual Deviance: 8381 AIC: 203.2
anova(x1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 12648.8
## shrub_density 1 2257.7 20 10391.1 0.02766 *
## aridity 1 1570.4 19 8820.6 0.06628 .
## shrub_density:aridity 1 439.8 18 8380.8 0.33110
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(x1), "\n")
## Linear Model AIC: 203.1718
quadratic_model_x1 <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian", data = final_data_2023)
quadratic_model_x1
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + shrub_density:aridity, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 68.7465 -5.2639 0.3959
## aridity I(aridity^2) shrub_density:aridity
## -17.6832 1.0897 0.2494
##
## Degrees of Freedom: 21 Total (i.e. Null); 16 Residual
## Null Deviance: 12650
## Residual Deviance: 7045 AIC: 203.4
anova(quadratic_model_x1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 12648.8
## shrub_density 1 2257.72 20 10391.1 0.02355 *
## I(shrub_density^2) 1 55.46 19 10335.6 0.72267
## aridity 1 1622.60 18 8713.0 0.05490 .
## I(aridity^2) 1 1189.72 17 7523.3 0.10022
## shrub_density:aridity 1 478.47 16 7044.8 0.29720
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_x1), "\n")
## Quadratic Model AIC: 203.3514
if (AIC(quadratic_model_x1) < AIC(x1)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
abund_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, animals)) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_dens_2023
x2 <- glm(richness ~ shrub_density * aridity, data = final_data_2023)
x2
##
## Call: glm(formula = richness ~ shrub_density * aridity, data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density aridity
## 1.79659 0.26232 0.28474
## shrub_density:aridity
## -0.00865
##
## Degrees of Freedom: 21 Total (i.e. Null); 18 Residual
## Null Deviance: 95.32
## Residual Deviance: 52.71 AIC: 91.66
anova(x2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 95.318
## shrub_density 1 25.4984 20 69.820 0.003169 **
## aridity 1 16.3249 19 53.495 0.018219 *
## shrub_density:aridity 1 0.7861 18 52.709 0.604382
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(x2), "\n")
## Linear Model AIC: 91.65558
quadratic_model_x2 <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
quadratic_model_x2
##
## Call: glm(formula = richness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 2.7345590 -0.0007208 0.0196188
## aridity I(aridity^2) shrub_density:aridity
## 0.0227789 0.0152190 -0.0050956
##
## Degrees of Freedom: 21 Total (i.e. Null); 16 Residual
## Null Deviance: 95.32
## Residual Deviance: 51.78 AIC: 95.26
anova(quadratic_model_x2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 95.318
## shrub_density 1 25.4984 20 69.820 0.00500 **
## I(shrub_density^2) 1 0.4031 19 69.417 0.72413
## aridity 1 17.3499 18 52.067 0.02059 *
## I(aridity^2) 1 0.0894 17 51.977 0.86800
## shrub_density:aridity 1 0.1998 16 51.778 0.80379
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_x2), "\n")
## Quadratic Model AIC: 95.26343
if (AIC(quadratic_model_x2) < AIC(x2)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
rich_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, richness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_dens_2023
x3 <- glm(Average_Evenness ~ shrub_density * aridity, data = final_data_2023)
x3
##
## Call: glm(formula = Average_Evenness ~ shrub_density * aridity, data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density aridity
## -0.0315882 0.0055497 0.0122207
## shrub_density:aridity
## -0.0001543
##
## Degrees of Freedom: 21 Total (i.e. Null); 18 Residual
## Null Deviance: 0.1608
## Residual Deviance: 0.1076 AIC: -44.6
anova(x3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.16080
## shrub_density 1 0.013715 20 0.14708 0.12993
## aridity 1 0.039182 19 0.10790 0.01048 *
## shrub_density:aridity 1 0.000250 18 0.10765 0.83796
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(x3), "\n")
## Linear Model AIC: -44.60481
quadratic_model_x3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
quadratic_model_x3
##
## Call: glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + shrub_density:aridity, data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## -7.236e-01 3.997e-03 -2.517e-05
## aridity I(aridity^2) shrub_density:aridity
## 2.134e-01 -1.181e-02 2.469e-04
##
## Degrees of Freedom: 21 Total (i.e. Null); 16 Residual
## Null Deviance: 0.1608
## Residual Deviance: 0.006624 AIC: -101.9
anova(quadratic_model_x3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.160797
## shrub_density 1 0.013715 20 0.147082 8.624e-09 ***
## I(shrub_density^2) 1 0.020159 19 0.126922 2.991e-12 ***
## aridity 1 0.023110 18 0.103812 7.927e-14 ***
## I(aridity^2) 1 0.096719 17 0.007093 < 2.2e-16 ***
## shrub_density:aridity 1 0.000469 16 0.006624 0.2872
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_model_x3), "\n")
## Quadratic Model AIC: -101.9453
if (AIC(quadratic_model_x1) < AIC(x1)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
even_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_dens_2023
plot2 <- abund_dens_2023/rich_dens_2023/even_dens_2023
plot2
### Temperature 2022 (Microsite/plot level)
t1 <- glm(animals ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2022)
t1
##
## Call: glm(formula = animals ~ shrub_density * mean_temp, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density mean_temp
## 641.3669 22.3664 -20.8447
## shrub_density:mean_temp
## -0.6221
##
## Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
## Null Deviance: 619900
## Residual Deviance: 389100 AIC: 310.8
anova(t1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 619927
## shrub_density 1 26978 22 592949 0.238985
## mean_temp 1 199846 21 393103 0.001351 **
## shrub_density:mean_temp 1 3970 20 389134 0.651496
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(t1), "\n")
## Linear Model AIC: 310.756
quadratic_temp_1 <- glm(animals ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2022)
quadratic_temp_1
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 4363.574 155.190 -6.801
## mean_temp I(mean_temp^2) shrub_density:mean_temp
## -305.916 5.317 -2.618
##
## Degrees of Freedom: 23 Total (i.e. Null); 18 Residual
## Null Deviance: 619900
## Residual Deviance: 301200 AIC: 308.6
anova(quadratic_temp_1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 619927
## shrub_density 1 26978 22 592949 0.2042037
## I(shrub_density^2) 1 2942 21 590007 0.6750263
## mean_temp 1 240180 20 349827 0.0001516 ***
## I(mean_temp^2) 1 2772 19 347055 0.6840117
## shrub_density:mean_temp 1 45820 18 301235 0.0979909 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_1), "\n")
## Quadratic Model AIC: 308.6113
if (AIC(quadratic_temp_1) < AIC(t1)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
abund_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Mean Temperature (°C)", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_temp_2022
t2 <- glm(richness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2022)
t2
##
## Call: glm(formula = richness ~ shrub_density * mean_temp, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density mean_temp
## 32.305966 -0.039811 -1.022221
## shrub_density:mean_temp
## 0.005855
##
## Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
## Null Deviance: 475
## Residual Deviance: 120.1 AIC: 116.8
anova(t2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 474.96
## shrub_density 1 6.43 22 468.53 0.3008
## mean_temp 1 348.04 21 120.49 2.703e-14 ***
## shrub_density:mean_temp 1 0.35 20 120.14 0.8088
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(t2), "\n")
## Linear Model AIC: 116.7632
quadratic_temp_2 <- glm(richness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2022)
quadratic_temp_2
##
## Call: glm(formula = richness ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 139.969166 -1.240892 0.104613
## mean_temp I(mean_temp^2) shrub_density:mean_temp
## -9.215104 0.152255 0.007732
##
## Degrees of Freedom: 23 Total (i.e. Null); 18 Residual
## Null Deviance: 475
## Residual Deviance: 74.51 AIC: 109.3
anova(quadratic_temp_2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 474.96
## shrub_density 1 6.431 22 468.53 0.21261
## I(shrub_density^2) 1 106.430 21 362.10 3.965e-07 ***
## mean_temp 1 262.720 20 99.38 1.631e-15 ***
## I(mean_temp^2) 1 24.467 19 74.91 0.01505 *
## shrub_density:mean_temp 1 0.400 18 74.51 0.75602
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_2), "\n")
## Quadratic Model AIC: 109.2983
if (AIC(quadratic_temp_2) < AIC(t2)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
rich_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Mean Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_temp_2022
t3 <- glm(Average_Evenness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2022)
t3
##
## Call: glm(formula = Average_Evenness ~ shrub_density * mean_temp, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density mean_temp
## 0.4144890 0.0186680 -0.0124555
## shrub_density:mean_temp
## -0.0005946
##
## Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
## Null Deviance: 0.1455
## Residual Deviance: 0.05211 AIC: -69.07
anova(t3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.145519
## shrub_density 1 0.007333 22 0.138186 0.09341 .
## mean_temp 1 0.082454 21 0.055732 1.847e-08 ***
## shrub_density:mean_temp 1 0.003627 20 0.052105 0.23807
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(t3), "\n")
## Linear Model AIC: -69.07199
quadratic_temp_3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2022)
quadratic_temp_3
##
## Call: glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian",
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 1.8333330 -0.0250452 0.0030448
## mean_temp I(mean_temp^2) shrub_density:mean_temp
## -0.1201345 0.0019981 -0.0002403
##
## Degrees of Freedom: 23 Total (i.e. Null); 18 Residual
## Null Deviance: 0.1455
## Residual Deviance: 0.02914 AIC: -79.02
anova(quadratic_temp_3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.145519
## shrub_density 1 0.007333 22 0.138186 0.03332 *
## I(shrub_density^2) 1 0.056383 21 0.081804 3.601e-09 ***
## mean_temp 1 0.049707 20 0.032096 3.004e-08 ***
## I(mean_temp^2) 1 0.002571 19 0.029526 0.20760
## shrub_density:mean_temp 1 0.000386 18 0.029140 0.62536
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_3), "\n")
## Quadratic Model AIC: -79.02002
if (AIC(quadratic_temp_3) < AIC(t3)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
even_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Mean Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2022
### Temperature 2023 (Microsite/plot level)
z1 <- glm(animals ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2023)
z1
##
## Call: glm(formula = animals ~ shrub_density * mean_temp, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density mean_temp
## 18.0908 12.3265 -0.1743
## shrub_density:mean_temp
## -0.4532
##
## Degrees of Freedom: 21 Total (i.e. Null); 18 Residual
## Null Deviance: 12650
## Residual Deviance: 8513 AIC: 203.5
anova(z1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 12648.8
## shrub_density 1 2257.72 20 10391.1 0.02890 *
## mean_temp 1 1313.77 19 9077.3 0.09558 .
## shrub_density:mean_temp 1 563.83 18 8513.5 0.27491
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(z1), "\n")
## Linear Model AIC: 203.5172
quadratic_temp_z1 <- glm(animals ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2023)
quadratic_temp_z1
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 1703.91419 4.63774 -0.02144
## mean_temp I(mean_temp^2) shrub_density:mean_temp
## -139.96488 2.88050 -0.14516
##
## Degrees of Freedom: 21 Total (i.e. Null); 16 Residual
## Null Deviance: 12650
## Residual Deviance: 5706 AIC: 198.7
anova(quadratic_temp_z1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 12648.8
## shrub_density 1 2257.72 20 10391.1 0.011869 *
## I(shrub_density^2) 1 55.46 19 10335.6 0.693345
## mean_temp 1 1517.38 18 8818.2 0.039146 *
## I(mean_temp^2) 1 3064.70 17 5753.5 0.003375 **
## shrub_density:mean_temp 1 47.11 16 5706.4 0.716289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_z1), "\n")
## Quadratic Model AIC: 198.716
if (AIC(quadratic_temp_z1) < AIC(z1)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
abund_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Mean Temperature (°C)", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_temp_2023
z2 <- glm(richness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2023)
z2
##
## Call: glm(formula = richness ~ shrub_density * mean_temp, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density mean_temp
## -3.72064 0.97821 0.33694
## shrub_density:mean_temp
## -0.03281
##
## Degrees of Freedom: 21 Total (i.e. Null); 18 Residual
## Null Deviance: 95.32
## Residual Deviance: 66 AIC: 96.6
anova(z2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 95.318
## shrub_density 1 25.4984 20 69.820 0.008364 **
## mean_temp 1 0.8608 19 68.959 0.628021
## shrub_density:mean_temp 1 2.9560 18 66.003 0.369258
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(z2), "\n")
## Linear Model AIC: 96.60373
quadratic_temp_z2 <- glm(richness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2023)
quadratic_temp_z2
##
## Call: glm(formula = richness ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 66.1097226 0.6326525 0.0007925
## mean_temp I(mean_temp^2) shrub_density:mean_temp
## -5.4527083 0.1192890 -0.0197340
##
## Degrees of Freedom: 21 Total (i.e. Null); 16 Residual
## Null Deviance: 95.32
## Residual Deviance: 61.32 AIC: 98.99
anova(quadratic_temp_z2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 95.318
## shrub_density 1 25.4984 20 69.820 0.0099 **
## I(shrub_density^2) 1 0.4031 19 69.417 0.7457
## mean_temp 1 0.6475 18 68.769 0.6811
## I(mean_temp^2) 1 6.5745 17 62.195 0.1903
## shrub_density:mean_temp 1 0.8706 16 61.324 0.6337
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_z2), "\n")
## Quadratic Model AIC: 98.98617
if (AIC(quadratic_temp_z2) < AIC(z2)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Linear model is better or equally good.
rich_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Mean Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_temp_2023
z3 <- glm(Average_Evenness ~ shrub_density * mean_temp, family = "gaussian", data = final_data_2023)
z3
##
## Call: glm(formula = Average_Evenness ~ shrub_density * mean_temp, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density mean_temp
## -0.879560 0.033490 0.039643
## shrub_density:mean_temp
## -0.001127
##
## Degrees of Freedom: 21 Total (i.e. Null); 18 Residual
## Null Deviance: 0.1608
## Residual Deviance: 0.04301 AIC: -64.79
anova(z3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.160797
## shrub_density 1 0.013715 20 0.147082 0.01659 *
## mean_temp 1 0.100583 19 0.046499 8.713e-11 ***
## shrub_density:mean_temp 1 0.003485 18 0.043014 0.22718
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(z3), "\n")
## Linear Model AIC: -64.78675
quadratic_temp_z3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian", data = final_data_2023)
quadratic_temp_z3
##
## Call: glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + shrub_density:mean_temp, family = "gaussian",
## data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2)
## 5.642e+00 4.705e-03 -1.426e-04
## mean_temp I(mean_temp^2) shrub_density:mean_temp
## -5.012e-01 1.114e-02 5.355e-05
##
## Degrees of Freedom: 21 Total (i.e. Null); 16 Residual
## Null Deviance: 0.1608
## Residual Deviance: 0.0005175 AIC: -158
anova(quadratic_temp_z3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.160797
## shrub_density 1 0.013715 20 0.147082 <2e-16 ***
## I(shrub_density^2) 1 0.020159 19 0.126922 <2e-16 ***
## mean_temp 1 0.085588 18 0.041335 <2e-16 ***
## I(mean_temp^2) 1 0.040811 17 0.000524 <2e-16 ***
## shrub_density:mean_temp 1 0.000006 16 0.000518 0.6562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(quadratic_temp_z3), "\n")
## Quadratic Model AIC: -158.0316
if (AIC(quadratic_temp_z3) < AIC(z3)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
even_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Mean Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2023
### Test colinearity 2022
your_model <- glm(animals ~ shrub_density + I(shrub_density^2) + mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity + I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
your_model
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity +
## I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) mean_temp
## 68.06110 0.63570 -0.05244 6.70852
## I(mean_temp^2) aridity I(aridity^2) mean_humidity
## -0.14576 -74.57567 7.80189 2.59059
## I(mean_humidity^2)
## -0.04495
##
## Degrees of Freedom: 23 Total (i.e. Null); 15 Residual
## Null Deviance: 4036
## Residual Deviance: 416.3 AIC: 562.1
summary(your_model)
##
## Call:
## glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity +
## I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 68.06110 30.34372 2.243 0.024897 *
## shrub_density 0.63570 0.02693 23.603 < 2e-16 ***
## I(shrub_density^2) -0.05244 0.00238 -22.037 < 2e-16 ***
## mean_temp 6.70852 1.96996 3.405 0.000661 ***
## I(mean_temp^2) -0.14576 0.04251 -3.429 0.000607 ***
## aridity -74.57567 21.89055 -3.407 0.000657 ***
## I(aridity^2) 7.80189 2.34090 3.333 0.000860 ***
## mean_humidity 2.59059 0.75099 3.450 0.000561 ***
## I(mean_humidity^2) -0.04495 0.01281 -3.509 0.000450 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4036.16 on 23 degrees of freedom
## Residual deviance: 416.33 on 15 degrees of freedom
## AIC: 562.14
##
## Number of Fisher Scoring iterations: 6
your_model
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## mean_temp + I(mean_temp^2) + aridity + I(aridity^2) + mean_humidity +
## I(mean_humidity^2), family = poisson(link = "log"), data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) mean_temp
## 68.06110 0.63570 -0.05244 6.70852
## I(mean_temp^2) aridity I(aridity^2) mean_humidity
## -0.14576 -74.57567 7.80189 2.59059
## I(mean_humidity^2)
## -0.04495
##
## Degrees of Freedom: 23 Total (i.e. Null); 15 Residual
## Null Deviance: 4036
## Residual Deviance: 416.3 AIC: 562.1
a <- anova(your_model, test = "Chisq")
a
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 4036.2
## shrub_density 1 205.26 22 3830.9 < 2.2e-16 ***
## I(shrub_density^2) 1 22.68 21 3808.2 1.91e-06 ***
## mean_temp 1 2700.17 20 1108.0 < 2.2e-16 ***
## I(mean_temp^2) 1 0.00 19 1108.0 0.9503221
## aridity 1 669.01 18 439.0 < 2.2e-16 ***
## I(aridity^2) 1 3.95 17 435.1 0.0468959 *
## mean_humidity 1 6.48 16 428.6 0.0109107 *
## I(mean_humidity^2) 1 12.27 15 416.3 0.0004594 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(car)
vif_results <- car::vif(your_model)
print(vif_results)
## shrub_density I(shrub_density^2) mean_temp I(mean_temp^2)
## 68.80459 72.66716 24431.18644 29383.10054
## aridity I(aridity^2) mean_humidity I(mean_humidity^2)
## 235991.70543 224724.98700 15974.68265 10607.89891
cor_matrix <- cor(final_data_2022[, c("shrub_density", "mean_temp", "aridity", "mean_humidity")])
print(cor_matrix)
## shrub_density mean_temp aridity mean_humidity
## shrub_density 1.00000000 0.02938781 0.02585627 0.04761062
## mean_temp 0.02938781 1.00000000 -0.85229523 -0.98132444
## aridity 0.02585627 -0.85229523 1.00000000 0.80272028
## mean_humidity 0.04761062 -0.98132444 0.80272028 1.00000000
library(lmerTest)
library(performance)
model <- lmer(animals ~ shrub_density + mean_temp + mean_humidity + (1|year), data = final_data)
check_collinearity(model)
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## shrub_density 1.11 [1.01, 2.77] 1.05 0.90 [0.36, 0.99]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## mean_temp 7.50 [4.80, 12.13] 2.74 0.13 [0.08, 0.21]
## mean_humidity 7.68 [4.90, 12.42] 2.77 0.13 [0.08, 0.20]
anova(model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## shrub_density 26141 26141 1 41.354 2.4335 0.12639
## mean_temp 61743 61743 1 41.634 5.7477 0.02108 *
## mean_humidity 12699 12699 1 39.750 1.1822 0.28347
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ranova(model)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## animals ~ shrub_density + mean_temp + mean_humidity + (1 | year)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 -267.52 547.03
## (1 | year) 5 -269.78 549.56 4.5275 1 0.03335 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### So use one model for everything!
### 2022 All together since correlation is low/moderate
a1 <- glm(animals ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2022)
a1
##
## Call: glm(formula = animals ~ shrub_density + aridity + mean_temp +
## mean_humidity, data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density aridity mean_temp mean_humidity
## 594.701 5.718 -87.742 -15.161 12.542
##
## Degrees of Freedom: 23 Total (i.e. Null); 19 Residual
## Null Deviance: 619900
## Residual Deviance: 339100 AIC: 309.5
anova(a1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 619927
## shrub_density 1 26978 22 592949 0.218915
## aridity 1 71496 21 521453 0.045349 *
## mean_temp 1 176488 20 344964 0.001664 **
## mean_humidity 1 5834 19 339131 0.567519
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a1), "\n")
## Linear Model AIC: 309.4551
#q_model_1 <- glmer(animals ~ shrub_density + aridity + mean_temp, data = final_data_2022)
#q_model_1
#anova(q_model_1, test = "Chisq")
#cat("Quadratic Model AIC:", AIC(q_model_1), "\n")
if (AIC(quadratic_model_1) < AIC(a1)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
a2 <- glm(richness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2022)
a2
##
## Call: glm(formula = richness ~ shrub_density + aridity + mean_temp +
## mean_humidity, data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density aridity mean_temp mean_humidity
## 5.01364 0.05630 -3.01765 -0.08542 0.73972
##
## Degrees of Freedom: 23 Total (i.e. Null); 19 Residual
## Null Deviance: 475
## Residual Deviance: 29.85 AIC: 85.35
anova(a2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 474.96
## shrub_density 1 6.431 22 468.53 0.0430701 *
## aridity 1 133.613 21 334.91 < 2.2e-16 ***
## mean_temp 1 284.765 20 50.15 < 2.2e-16 ***
## mean_humidity 1 20.295 19 29.85 0.0003258 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a2), "\n")
## Linear Model AIC: 85.34806
q_model_2 <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2022)
q_model_2
##
## Call: glm(formula = richness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity +
## I(mean_humidity^2), data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) aridity
## 107.7409 -0.8030 0.0688 99.3440
## I(aridity^2) mean_temp I(mean_temp^2) mean_humidity
## -10.8127 -37.2859 0.8090 7.0335
## I(mean_humidity^2)
## -0.1173
##
## Degrees of Freedom: 23 Total (i.e. Null); 15 Residual
## Null Deviance: 475
## Residual Deviance: 15.68 AIC: 77.89
anova(q_model_2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 474.96
## shrub_density 1 6.431 22 468.53 0.013127 *
## I(shrub_density^2) 1 106.430 21 362.10 < 2.2e-16 ***
## aridity 1 101.019 20 261.08 < 2.2e-16 ***
## I(aridity^2) 1 232.882 19 28.20 < 2.2e-16 ***
## mean_temp 1 0.134 18 28.06 0.720348
## I(mean_temp^2) 1 10.796 17 17.27 0.001311 **
## mean_humidity 1 0.371 16 16.89 0.551123
## I(mean_humidity^2) 1 1.215 15 15.68 0.281016
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_2), "\n")
## Quadratic Model AIC: 77.89283
if (AIC(q_model_2) < AIC(a2)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
a3 <- glm(Average_Evenness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2022)
a3
##
## Call: glm(formula = Average_Evenness ~ shrub_density + aridity + mean_temp +
## mean_humidity, data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density aridity mean_temp mean_humidity
## -0.983958 0.001086 -0.033891 0.026940 0.024481
##
## Degrees of Freedom: 23 Total (i.e. Null); 19 Residual
## Null Deviance: 0.1455
## Residual Deviance: 0.01438 AIC: -97.98
anova(a3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.145519
## shrub_density 1 0.007333 22 0.138186 0.001851 **
## aridity 1 0.029971 21 0.108215 3.098e-10 ***
## mean_temp 1 0.071611 20 0.036605 < 2.2e-16 ***
## mean_humidity 1 0.022229 19 0.014376 5.950e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a3), "\n")
## Linear Model AIC: -97.97747
q_model_3 <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2022)
q_model_3
##
## Call: glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity +
## I(mean_humidity^2), data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) aridity
## 2.766661 -0.026205 0.002274 -1.027071
## I(aridity^2) mean_temp I(mean_temp^2) mean_humidity
## 0.107674 -0.241692 0.005770 0.146552
## I(mean_humidity^2)
## -0.002244
##
## Degrees of Freedom: 23 Total (i.e. Null); 15 Residual
## Null Deviance: 0.1455
## Residual Deviance: 0.005106 AIC: -114.8
anova(q_model_3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.145519
## shrub_density 1 0.007333 22 0.138186 3.460e-06 ***
## I(shrub_density^2) 1 0.056383 21 0.081804 < 2.2e-16 ***
## aridity 1 0.018828 20 0.062975 1.026e-13 ***
## I(aridity^2) 1 0.049621 19 0.013354 < 2.2e-16 ***
## mean_temp 1 0.001140 18 0.012214 0.067244 .
## I(mean_temp^2) 1 0.003081 17 0.009133 0.002626 **
## mean_humidity 1 0.003583 16 0.005550 0.001176 **
## I(mean_humidity^2) 1 0.000444 15 0.005106 0.253210
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_3), "\n")
## Quadratic Model AIC: -114.822
if (AIC(q_model_3) < AIC(a3)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
b1 <- glm(animals ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2023)
b1
##
## Call: glm(formula = animals ~ shrub_density + aridity + mean_temp +
## mean_humidity, data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density aridity mean_temp mean_humidity
## 194.4432 1.3948 3.4741 -8.0000 -0.6446
##
## Degrees of Freedom: 21 Total (i.e. Null); 17 Residual
## Null Deviance: 12650
## Residual Deviance: 5738 AIC: 196.8
anova(b1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 12648.8
## shrub_density 1 2257.72 20 10391.1 0.009702 **
## aridity 1 1570.44 19 8820.6 0.031007 *
## mean_temp 1 3018.31 18 5802.3 0.002787 **
## mean_humidity 1 64.11 17 5738.2 0.662961
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(a1), "\n")
## Linear Model AIC: 309.4551
q_model_1b <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2023)
q_model_1b
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity +
## I(mean_humidity^2), data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) aridity
## -486.3157 -10.0215 0.7863 224.2800
## I(aridity^2) mean_temp I(mean_temp^2) mean_humidity
## -12.2965 -73.3075 1.4164 30.5245
## I(mean_humidity^2)
## -0.3424
##
## Degrees of Freedom: 21 Total (i.e. Null); 13 Residual
## Null Deviance: 12650
## Residual Deviance: 2388 AIC: 185.5
anova(q_model_1b, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 12648.8
## shrub_density 1 2257.72 20 10391.1 0.0004549 ***
## I(shrub_density^2) 1 55.46 19 10335.6 0.5826786
## aridity 1 1622.60 18 8713.0 0.0029563 **
## I(aridity^2) 1 1189.72 17 7523.3 0.0109254 *
## mean_temp 1 2572.77 16 4950.5 0.0001821 ***
## I(mean_temp^2) 1 1451.04 15 3499.5 0.0049430 **
## mean_humidity 1 529.64 14 2969.8 0.0894861 .
## I(mean_humidity^2) 1 582.08 13 2387.7 0.0750425 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_1b), "\n")
## Quadratic Model AIC: 185.5486
if (AIC(q_model_1b) < AIC(b1)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
b2 <- glm(richness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2023)
b2
##
## Call: glm(formula = richness ~ shrub_density + aridity + mean_temp +
## mean_humidity, data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density aridity mean_temp mean_humidity
## -13.0027 0.1575 0.2489 0.4385 0.1494
##
## Degrees of Freedom: 21 Total (i.e. Null); 17 Residual
## Null Deviance: 95.32
## Residual Deviance: 49.66 AIC: 92.35
anova(b2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 95.318
## shrub_density 1 25.4984 20 69.820 0.003133 **
## aridity 1 16.3249 19 53.495 0.018084 *
## mean_temp 1 0.3848 18 53.110 0.716651
## mean_humidity 1 3.4457 17 49.664 0.277465
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(b2), "\n")
## Linear Model AIC: 92.34667
q_model_2b <- glm(richness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2023)
q_model_2b
##
## Call: glm(formula = richness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity +
## I(mean_humidity^2), data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) aridity
## -172.23925 -0.98555 0.08180 24.28498
## I(aridity^2) mean_temp I(mean_temp^2) mean_humidity
## -1.33580 2.07678 -0.04053 2.86793
## I(mean_humidity^2)
## -0.02884
##
## Degrees of Freedom: 21 Total (i.e. Null); 13 Residual
## Null Deviance: 95.32
## Residual Deviance: 26.03 AIC: 86.13
anova(q_model_2b, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 95.318
## shrub_density 1 25.4984 20 69.820 0.0003585 ***
## I(shrub_density^2) 1 0.4031 19 69.417 0.6536249
## aridity 1 17.3499 18 52.067 0.0032410 **
## I(aridity^2) 1 0.0894 17 51.977 0.8326496
## mean_temp 1 0.3069 16 51.670 0.6953844
## I(mean_temp^2) 1 8.5804 15 43.090 0.0384262 *
## mean_humidity 1 12.9355 14 30.155 0.0110236 *
## I(mean_humidity^2) 1 4.1294 13 26.025 0.1509434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_2b), "\n")
## Quadratic Model AIC: 86.12974
if (AIC(q_model_2b) < AIC(b2)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
b3 <- glm(Average_Evenness ~ shrub_density + aridity + mean_temp + mean_humidity, data = final_data_2023)
b3
##
## Call: glm(formula = Average_Evenness ~ shrub_density + aridity + mean_temp +
## mean_humidity, data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density aridity mean_temp mean_humidity
## 0.612360 0.008000 0.004847 -0.009961 -0.011192
##
## Degrees of Freedom: 21 Total (i.e. Null); 17 Residual
## Null Deviance: 0.1608
## Residual Deviance: 0.01975 AIC: -79.91
anova(b3, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.160797
## shrub_density 1 0.013715 20 0.147082 0.0005914 ***
## aridity 1 0.039182 19 0.107900 6.369e-09 ***
## mean_temp 1 0.068815 18 0.039085 1.410e-14 ***
## mean_humidity 1 0.019331 17 0.019754 4.530e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Linear Model AIC:", AIC(b3), "\n")
## Linear Model AIC: -79.90611
q_model_3b <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity + I(mean_humidity^2), data = final_data_2023)
q_model_3b
##
## Call: glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + mean_temp + I(mean_temp^2) + mean_humidity +
## I(mean_humidity^2), data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) aridity
## 6.433e+00 5.868e-03 -1.715e-04 -3.056e-02
## I(aridity^2) mean_temp I(mean_temp^2) mean_humidity
## 1.657e-03 -5.691e-01 1.279e-02 -1.462e-03
## I(mean_humidity^2)
## 4.833e-05
##
## Degrees of Freedom: 21 Total (i.e. Null); 13 Residual
## Null Deviance: 0.1608
## Residual Deviance: 0.0002097 AIC: -171.9
anova(q_model_3b, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.160797
## shrub_density 1 0.013715 20 0.147082 < 2.2e-16 ***
## I(shrub_density^2) 1 0.020159 19 0.126922 < 2.2e-16 ***
## aridity 1 0.023110 18 0.103812 < 2.2e-16 ***
## I(aridity^2) 1 0.096719 17 0.007093 < 2.2e-16 ***
## mean_temp 1 0.001253 16 0.005840 < 2.2e-16 ***
## I(mean_temp^2) 1 0.005458 15 0.000381 < 2.2e-16 ***
## mean_humidity 1 0.000160 14 0.000221 0.001638 **
## I(mean_humidity^2) 1 0.000012 13 0.000210 0.396490
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("Quadratic Model AIC:", AIC(q_model_3b), "\n")
## Quadratic Model AIC: -171.9054
if (AIC(q_model_3b) < AIC(b3)) {
cat("Quadratic model is better.\n")
} else {
cat("Linear model is better or equally good.\n")
}
## Quadratic model is better.
library(sjPlot)
q_model_1 <- lmer(Average_Evenness ~ shrub_density + aridity + mean_temp + (1|year), data = final_data) ### Generalize linear mixed effect models
q_model_1
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: Average_Evenness ~ shrub_density + aridity + mean_temp + (1 |
## year)
## Data: final_data
## REML criterion at convergence: -74.5199
## Random effects:
## Groups Name Std.Dev.
## year (Intercept) 0.05073
## Residual 0.07346
## Number of obs: 46, groups: year, 2
## Fixed Effects:
## (Intercept) shrub_density aridity mean_temp
## 0.120540 0.003439 0.011743 -0.004518
anova(q_model_1, test = "Chisq")
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## shrub_density 0.017489 0.017489 1 41.000 3.2410 0.079177 .
## aridity 0.048741 0.048741 1 36.939 9.0326 0.004745 **
## mean_temp 0.009448 0.009448 1 41.643 1.7509 0.192988
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("LMER Model AIC:", AIC(q_model_1), "\n")
## LMER Model AIC: -62.51993
tab_model(q_model_1)
|
Â
|
Average_Evenness
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
0.12
|
-0.08 – 0.32
|
0.228
|
|
shrub density
|
0.00
|
-0.00 – 0.01
|
0.079
|
|
aridity
|
0.01
|
0.00 – 0.02
|
0.005
|
|
mean temp
|
-0.00
|
-0.01 – 0.00
|
0.193
|
|
Random Effects
|
|
σ2
|
0.01
|
|
Ï„00 year
|
0.00
|
|
ICC
|
0.32
|
|
N year
|
2
|
|
Observations
|
46
|
|
Marginal R2 / Conditional R2
|
0.276 / 0.510
|
library(ggplot2)
even_temp <- ggplot(final_data, aes(aridity, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ x +I(x^2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "aridity", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp
cor.test(final_data_2022$aridity, final_data_2022$mean_humidity,method = c("pearson"))
##
## Pearson's product-moment correlation
##
## data: final_data_2022$aridity and final_data_2022$mean_humidity
## t = 6.3135, df = 22, p-value = 2.36e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5905535 0.9110920
## sample estimates:
## cor
## 0.8027203
q_model_2 <- lmer(animals ~ shrub_density + aridity + mean_temp + (1|year), data = final_data) ### Generalize linear mixed effect models
q_model_2
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: animals ~ shrub_density + aridity + mean_temp + (1 | year)
## Data: final_data
## REML criterion at convergence: 535.3932
## Random effects:
## Groups Name Std.Dev.
## year (Intercept) 123.4
## Residual 103.4
## Number of obs: 46, groups: year, 2
## Fixed Effects:
## (Intercept) shrub_density aridity mean_temp
## 489.791 3.422 4.482 -18.519
anova(q_model_2, test = "Chisq")
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## shrub_density 17315 17315 1 41.000 1.621 0.2101201
## aridity 6761 6761 1 41.958 0.633 0.4307392
## mean_temp 157818 157818 1 41.275 14.775 0.0004111 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("LMER Model AIC:", AIC(q_model_2), "\n")
## LMER Model AIC: 547.3932
tab_model(q_model_2)
|
Â
|
animals
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
489.79
|
175.18 – 804.41
|
0.003
|
|
shrub density
|
3.42
|
-2.01 – 8.85
|
0.210
|
|
aridity
|
4.48
|
-6.90 – 15.87
|
0.431
|
|
mean temp
|
-18.52
|
-28.26 – -8.78
|
<0.001
|
|
Random Effects
|
|
σ2
|
10681.71
|
|
Ï„00 year
|
15221.82
|
|
ICC
|
0.59
|
|
N year
|
2
|
|
Observations
|
46
|
|
Marginal R2 / Conditional R2
|
0.169 / 0.657
|
q_model_3 <- lmer(richness ~ shrub_density + aridity + mean_temp + mean_humidity + (1|year), data = final_data) ### Generalize linear mixed effect models
q_model_3
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: richness ~ shrub_density + aridity + mean_temp + mean_humidity +
## (1 | year)
## Data: final_data
## REML criterion at convergence: 222.1085
## Random effects:
## Groups Name Std.Dev.
## year (Intercept) 2.725
## Residual 2.449
## Number of obs: 46, groups: year, 2
## Fixed Effects:
## (Intercept) shrub_density aridity mean_temp mean_humidity
## 27.11280 0.12703 0.32017 -0.89528 -0.06545
anova(q_model_3, test = "Chisq")
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## shrub_density 21.243 21.243 1 40.459 3.5429 0.067006 .
## aridity 32.049 32.049 1 38.813 5.3452 0.026175 *
## mean_temp 45.401 45.401 1 40.386 7.5720 0.008835 **
## mean_humidity 1.494 1.494 1 38.798 0.2491 0.620509
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("LMER Model AIC:", AIC(q_model_3), "\n")
## LMER Model AIC: 236.1085
tab_model(q_model_3)
|
Â
|
richness
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
27.11
|
3.21 – 51.01
|
0.027
|
|
shrub density
|
0.13
|
-0.01 – 0.26
|
0.067
|
|
aridity
|
0.32
|
0.04 – 0.60
|
0.026
|
|
mean temp
|
-0.90
|
-1.55 – -0.24
|
0.009
|
|
mean humidity
|
-0.07
|
-0.33 – 0.20
|
0.620
|
|
Random Effects
|
|
σ2
|
6.00
|
|
Ï„00 year
|
7.43
|
|
ICC
|
0.55
|
|
N year
|
2
|
|
Observations
|
46
|
|
Marginal R2 / Conditional R2
|
0.414 / 0.738
|
library(emmeans)
arid <- glm(aridity~site * year, data = final_data)
summary(arid)
##
## Call:
## glm(formula = aridity ~ site * year, data = final_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.661e+04 2.570e-08 -6.463e+11 <2e-16 ***
## siteCuyama 1.614e+04 3.502e-08 4.608e+11 <2e-16 ***
## siteTecopa -1.536e+02 3.502e-08 -4.387e+09 <2e-16 ***
## year 8.215e+00 1.271e-11 6.466e+11 <2e-16 ***
## siteCuyama:year -7.980e+00 1.731e-11 -4.609e+11 <2e-16 ***
## siteTecopa:year 7.499e-02 1.731e-11 4.331e+09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 5.534469e-22)
##
## Null deviance: 6.4447e+02 on 45 degrees of freedom
## Residual deviance: 2.2138e-20 on 40 degrees of freedom
## AIC: -2113.4
##
## Number of Fisher Scoring iterations: 1
anova(arid, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: aridity
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 45 644.47
## site 2 137.97 43 506.50 < 2.2e-16 ***
## year 1 339.23 42 167.28 < 2.2e-16 ***
## site:year 2 167.28 40 0.00 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(arid, pairwise~site|year)
## $emmeans
## year = 2022:
## site emmean SE df lower.CL upper.CL
## Carrizo 5.13 8.32e-12 40 5.13 5.13
## Cuyama 4.35 8.32e-12 40 4.35 4.35
## Tecopa 3.12 8.32e-12 40 3.12 3.12
##
## year = 2023:
## site emmean SE df lower.CL upper.CL
## Carrizo 13.34 9.60e-12 40 13.34 13.34
## Cuyama 4.59 8.32e-12 40 4.59 4.59
## Tecopa 11.41 8.32e-12 40 11.41 11.41
##
## Confidence level used: 0.95
##
## $contrasts
## year = 2022:
## contrast estimate SE df t.ratio p.value
## Carrizo - Cuyama 0.778 1.18e-11 40 66120078152.000 <.0001
## Carrizo - Tecopa 2.010 1.18e-11 40 170900140809.000 <.0001
## Cuyama - Tecopa 1.232 1.18e-11 40 104780062301.000 <.0001
##
## year = 2023:
## contrast estimate SE df t.ratio p.value
## Carrizo - Cuyama 8.758 1.27e-11 40 689337834957.000 <.0001
## Carrizo - Tecopa 1.935 1.27e-11 40 152320861822.000 <.0001
## Cuyama - Tecopa -6.823 1.18e-11 40 -580044624656.000 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates
###
temp <- glm(mean_temp~site * year, data = final_data)
summary(temp)
##
## Call:
## glm(formula = mean_temp ~ site * year, data = final_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3422.1290 933.9302 3.664 0.000720 ***
## siteCuyama -2454.3538 1272.7531 -1.928 0.060927 .
## siteTecopa 6536.5176 1272.7531 5.136 7.68e-06 ***
## year -1.6806 0.4618 -3.639 0.000774 ***
## siteCuyama:year 1.2132 0.6293 1.928 0.060996 .
## siteTecopa:year -3.2291 0.6293 -5.131 7.79e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.7311318)
##
## Null deviance: 516.406 on 45 degrees of freedom
## Residual deviance: 29.245 on 40 degrees of freedom
## AIC: 123.71
##
## Number of Fisher Scoring iterations: 2
anova(temp, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: mean_temp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 45 516.41
## site 2 380.19 43 136.22 < 2.2e-16 ***
## year 1 65.07 42 71.15 < 2.2e-16 ***
## site:year 2 41.91 40 29.25 3.58e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(temp, pairwise~site|year)
## $emmeans
## year = 2022:
## site emmean SE df lower.CL upper.CL
## Carrizo 24.0 0.302 40 23.4 24.7
## Cuyama 22.8 0.302 40 22.2 23.4
## Tecopa 31.4 0.302 40 30.8 32.0
##
## year = 2023:
## site emmean SE df lower.CL upper.CL
## Carrizo 22.4 0.349 40 21.7 23.1
## Cuyama 22.3 0.302 40 21.7 22.9
## Tecopa 26.5 0.302 40 25.9 27.1
##
## Confidence level used: 0.95
##
## $contrasts
## year = 2022:
## contrast estimate SE df t.ratio p.value
## Carrizo - Cuyama 1.2539 0.428 40 2.933 0.0150
## Carrizo - Tecopa -7.3265 0.428 40 -17.137 <.0001
## Cuyama - Tecopa -8.5803 0.428 40 -20.069 <.0001
##
## year = 2023:
## contrast estimate SE df t.ratio p.value
## Carrizo - Cuyama 0.0407 0.462 40 0.088 0.9957
## Carrizo - Tecopa -4.0974 0.462 40 -8.873 <.0001
## Cuyama - Tecopa -4.1380 0.428 40 -9.679 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates
b1 <- glm(Average_Evenness ~ shrub_density + aridity + mean_temp, data = final_data_2022)
b1
##
## Call: glm(formula = Average_Evenness ~ shrub_density + aridity + mean_temp,
## data = final_data_2022)
##
## Coefficients:
## (Intercept) shrub_density aridity mean_temp
## 1.074552 0.003864 -0.065522 -0.027385
##
## Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
## Null Deviance: 0.1455
## Residual Deviance: 0.0366 AIC: -77.55
anova(b1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.145519
## shrub_density 1 0.007333 22 0.138186 0.04533 *
## aridity 1 0.029971 21 0.108215 5.195e-05 ***
## mean_temp 1 0.071611 20 0.036605 3.972e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Use second degree function for paper!
q <- glm(Average_Evenness ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
q
##
## Call: glm(formula = Average_Evenness ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) aridity
## 6.7398056 0.0063338 -0.0001901 -0.0516903
## I(aridity^2) mean_temp I(mean_temp^2)
## 0.0028418 -0.5858324 0.0130916
##
## Degrees of Freedom: 21 Total (i.e. Null); 15 Residual
## Null Deviance: 0.1608
## Residual Deviance: 0.0003813 AIC: -162.8
anova(q, test="Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.160797
## shrub_density 1 0.013715 20 0.147082 < 2.2e-16 ***
## I(shrub_density^2) 1 0.020159 19 0.126922 < 2.2e-16 ***
## aridity 1 0.023110 18 0.103812 < 2.2e-16 ***
## I(aridity^2) 1 0.096719 17 0.007093 < 2.2e-16 ***
## mean_temp 1 0.001253 16 0.005840 2.191e-12 ***
## I(mean_temp^2) 1 0.005458 15 0.000381 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
corelation_fig <- ggplot(final_data, aes(x = mean_temp, y = mean_humidity)) +
geom_point() + labs(fill = "", x = "Temperature (°C)", y = "Humidity (%)") +theme_classic()+theme(aspect.ratio = 1) + geom_smooth(method = lm, se = F) + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7)
cor <- corelation_fig + stat_cor(method = "pearson")
cor
### Inflection stuff
library(inflection)
q <- glm(animals ~ shrub_density + I(shrub_density^2) + aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
q
##
## Call: glm(formula = animals ~ shrub_density + I(shrub_density^2) +
## aridity + I(aridity^2) + mean_temp + I(mean_temp^2), data = final_data_2023)
##
## Coefficients:
## (Intercept) shrub_density I(shrub_density^2) aridity
## 3942.1978 -2.6344 0.2613 -114.4117
## I(aridity^2) mean_temp I(mean_temp^2)
## 6.4616 -309.4103 6.7500
##
## Degrees of Freedom: 21 Total (i.e. Null); 15 Residual
## Null Deviance: 12650
## Residual Deviance: 3499 AIC: 190
anova(q, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 12648.8
## shrub_density 1 2257.72 20 10391.1 0.0018654 **
## I(shrub_density^2) 1 55.46 19 10335.6 0.6258710
## aridity 1 1622.60 18 8713.0 0.0083581 **
## I(aridity^2) 1 1189.72 17 7523.3 0.0239314 *
## mean_temp 1 2572.77 16 4950.5 0.0008975 ***
## I(mean_temp^2) 1 1451.04 15 3499.5 0.0126336 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#abline(q, col = "red")
#inflection_point <- find_inflection(final_data$animals, predict(q))
#cat("Inflection Point:", inflection_point, "\n")
#inflection_point <- find_inflection(your_data$abundance, predict(model))
#cat("Inflection Point:", inflection_point, "\n")
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
### 2022 Density
library(inflection)
f=function(x){(1/8*x+1/2*x^2)/(1+1/8*x+1/2*x^2)}
x = final_data_2022$shrub_density
y = final_data_2022$animals
cc=check_curve(x,y)
cc
## $ctype
## [1] "concave"
##
## $index
## [1] 1
ese(x,y,cc$index)
## j1 j2 chi
## ESE 13 8 NaN
ipbese=bese(x,y,cc$index)
ipbese$iplast
## [1] NaN
inflection_point_dens_2022 <- final_data_2022[which.min(final_data_2022$y), ]
Fig1_a <- abund_dens_2022 + geom_vline(xintercept = inflection_point_dens_2022$shrub_density, linetype = 2)
Fig1_a
rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
theme(axis.title.x = element_blank()) + labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
# Assuming your data frame is named final_data_2022
# Assuming your data has columns 'shrub_density', 'animals', 'richness', 'evenness'
# Load necessary libraries
# Install and load the bcp package
#install.packages("changepoint")
library(changepoint)
# Function to find inflection point using bcp package
find_inflection_point_cpt <- function(x, y) {
data_ts <- data.frame(y = y)
cpt_result <- cpt.var(data_ts$y, method = "BinSeg")
change_point <- cpts(cpt_result)[1]
return(change_point)
}
# Fit models and find inflection points using bcp
inflection_animals <- find_inflection_point_cpt(final_data_2022$shrub_density, final_data_2022$animals)
inflection_richness <- find_inflection_point_cpt(final_data_2022$shrub_density, final_data_2022$richness)
inflection_evenness <- find_inflection_point_cpt(final_data_2022$shrub_density, final_data_2022$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", inflection_animals, "\n")
## Inflection Point for Animals: 12
cat("Inflection Point for Richness:", inflection_richness, "\n")
## Inflection Point for Richness: 8
cat("Inflection Point for Evenness:", inflection_evenness, "\n")
## Inflection Point for Evenness: NA
# Plot individual graphs
# ... (ggplot code and plotting)
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per20m Radius", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_a <- abund_dens_2022 + geom_vline(xintercept = inflection_animals, linetype = 2)
Fig1_a
rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_b <- rich_dens_2022 + geom_vline(xintercept = inflection_richness, linetype = 2)
Fig1_b
even_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") +
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_c <- even_dens_2022 + geom_vline(xintercept = inflection_evenness, linetype = 2)
Fig1_c
### Now for 2023
find_inflection_point_cpt <- function(x, y) {
data_ts <- data.frame(y = y)
cpt_result <- cpt.var(data_ts$y, method = "BinSeg")
change_point <- cpts(cpt_result)
return(change_point)
}
# Standardize variables
final_data_2022$standardized_aridity <- scale(final_data_2022$aridity)
# Adjust your model accordingly
inflection_animals <- find_inflection_point_cpt(final_data_2022$standardized_aridity, final_data_2022$animals)
inflection_richness <- find_inflection_point_cpt(final_data_2022$standardized_aridity, final_data_2022$richness)
inflection_evenness <- find_inflection_point_cpt(final_data_2022$standardized_aridity, final_data_2022$Average_Evenness)
# Fit models and find inflection points using changepoint
inflection_animals <- find_inflection_point_cpt(final_data_2022$aridity, final_data_2022$animals)
inflection_richness <- find_inflection_point_cpt(final_data_2022$aridity, final_data_2022$richness)
inflection_evenness <- find_inflection_point_cpt(final_data_2022$aridity, final_data_2022$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", ifelse(length(inflection_animals) > 0, inflection_animals, "No clear inflection point"), "\n")
## Inflection Point for Animals: 12
cat("Inflection Point for Richness:", ifelse(length(inflection_richness) > 0, inflection_richness, "No clear inflection point"), "\n")
## Inflection Point for Richness: 8
cat("Inflection Point for Evenness:", ifelse(length(inflection_evenness) > 0, inflection_evenness, "No clear inflection point"), "\n")
## Inflection Point for Evenness: No clear inflection point
# Function to find inflection point using quadratic model
find_inflection_point_quadratic <- function(x, y) {
quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
return(vertex)
}
# Find inflection points using quadratic model
den_inflection_animals <- find_inflection_point_quadratic(final_data_2022$shrub_density, final_data_2022$animals)
den_inflection_richness <- find_inflection_point_quadratic(final_data_2022$shrub_density, final_data_2022$richness)
den_inflection_evenness <- find_inflection_point_quadratic(final_data_2022$shrub_density, final_data_2022$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", den_inflection_animals, "\n")
## Inflection Point for Animals: 8.518166
cat("Inflection Point for Richness:", den_inflection_richness, "\n")
## Inflection Point for Richness: 5.595659
cat("Inflection Point for Evenness:", den_inflection_evenness, "\n")
## Inflection Point for Evenness: 5.493175
abund_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_a <- abund_dens_2022 + geom_vline(xintercept = den_inflection_animals, linetype = 2)
Fig1_a
rich_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_b <- rich_dens_2022 + geom_vline(xintercept = den_inflection_richness, linetype = 2)
Fig1_b
even_dens_2022 <- ggplot(final_data_2022, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_c <- even_dens_2022 + geom_vline(xintercept = den_inflection_evenness, linetype = 2)
Fig1_c
### Aridity
find_inflection_point_quadratic <- function(x, y) {
quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
return(vertex)
}
# Find inflection points using quadratic model
arid_inflection_animals <- find_inflection_point_quadratic(final_data_2022$aridity, final_data_2022$animals)
arid_inflection_richness <- find_inflection_point_quadratic(final_data_2022$aridity, final_data_2022$richness)
arid_inflection_evenness <- find_inflection_point_quadratic(final_data_2022$aridity, final_data_2022$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", inflection_animals, "\n")
## Inflection Point for Animals: 12
cat("Inflection Point for Richness:", inflection_richness, "\n")
## Inflection Point for Richness: 8
cat("Inflection Point for Evenness:", inflection_evenness, "\n")
## Inflection Point for Evenness:
abund_arid_2022 <- ggplot(final_data_2022, aes(aridity, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Aridity", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig2_a <- abund_arid_2022 + geom_vline(xintercept = arid_inflection_animals, linetype = 2)
Fig2_a
rich_arid_2022 <- ggplot(final_data_2022, aes(aridity, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Aridity", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig2_b <- rich_arid_2022 + geom_vline(xintercept = arid_inflection_richness, linetype = 2)
Fig2_b
even_arid_2022 <- ggplot(final_data_2022, aes(aridity, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Aridity", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig2_c <- even_arid_2022 + geom_vline(xintercept = arid_inflection_evenness, linetype = 2)
Fig2_c
# Temp
find_inflection_point_quadratic <- function(x, y) {
quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
return(vertex)
}
# Find inflection points using quadratic model
temp_inflection_animals_2022 <- find_inflection_point_quadratic(final_data_2022$mean_temp, final_data_2022$animals)
temp_inflection_richness_2022 <- find_inflection_point_quadratic(final_data_2022$mean_temp, final_data_2022$richness)
temp_inflection_evenness_2022 <- find_inflection_point_quadratic(final_data_2022$mean_temp, final_data_2022$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", temp_inflection_animals_2022, "\n")
## Inflection Point for Animals: 3.261247
cat("Inflection Point for Richness:", temp_inflection_richness_2022, "\n")
## Inflection Point for Richness: 31.81689
cat("Inflection Point for Evenness:", temp_inflection_evenness_2022, "\n")
## Inflection Point for Evenness: 46.82795
### Temperature might be linear not second degree?
abund_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, animals)) +
geom_smooth(method = "lm", se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Temperature (°C)", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
abund_temp_2022
Fig3_a <- abund_temp_2022
rich_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, richness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "B") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
rich_temp_2022
Fig3_b <- rich_temp_2022
even_temp_2022 <- ggplot(final_data_2022, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "C") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2022
Fig3_c <- even_temp_2022
find_inflection_point_quadratic <- function(x, y) {
quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
return(vertex)
}
# Find inflection points using quadratic model
den_inflection_animals_2023 <- find_inflection_point_quadratic(final_data_2023$shrub_density, final_data_2023$animals)
den_inflection_richness_2023 <- find_inflection_point_quadratic(final_data_2023$shrub_density, final_data_2023$richness)
den_inflection_evenness_2023 <- find_inflection_point_quadratic(final_data_2023$shrub_density, final_data_2023$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", den_inflection_animals_2023, "\n")
## Inflection Point for Animals: 11.78866
cat("Inflection Point for Richness:", den_inflection_richness_2023, "\n")
## Inflection Point for Richness: 13.25782
cat("Inflection Point for Evenness:", den_inflection_evenness_2023, "\n")
## Inflection Point for Evenness: 6.598329
abund_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "D") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_d <- abund_dens_2023 + geom_vline(xintercept = den_inflection_animals_2023, linetype = 2)
Fig1_d
rich_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) +
labs(tag = "E") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_e <- rich_dens_2023 + geom_vline(xintercept = den_inflection_richness_2023, linetype = 2)
Fig1_e
even_dens_2023 <- ggplot(final_data_2023, aes(shrub_density, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Shrub Density per 20m Radius", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+
labs(tag = "F") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig1_f <- even_dens_2023 + geom_vline(xintercept = den_inflection_evenness_2023, linetype = 2)
Fig1_f
### Aridity
find_inflection_point_quadratic <- function(x, y) {
quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
return(vertex)
}
# Find inflection points using quadratic model
arid_inflection_animals_2023 <- find_inflection_point_quadratic(final_data_2023$aridity, final_data_2023$animals)
arid_inflection_richness_2023 <- find_inflection_point_quadratic(final_data_2023$aridity, final_data_2023$richness)
arid_inflection_evenness_2023 <- find_inflection_point_quadratic(final_data_2023$aridity, final_data_2023$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", arid_inflection_animals_2023, "\n")
## Inflection Point for Animals: 7.870704
cat("Inflection Point for Richness:", arid_inflection_richness_2023, "\n")
## Inflection Point for Richness: 5.866049
cat("Inflection Point for Evenness:", arid_inflection_evenness_2023, "\n")
## Inflection Point for Evenness: 9.157136
abund_arid_2022 <- ggplot(final_data_2023, aes(aridity, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Aridity", y = "Animal Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig2_d <- abund_arid_2022 + geom_vline(xintercept = arid_inflection_animals_2023, linetype = 2)
Fig2_d
rich_arid_2022 <- ggplot(final_data_2023, aes(aridity, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Aridity", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig2_e <- rich_arid_2022 + geom_vline(xintercept = arid_inflection_richness_2023, linetype = 2)
Fig2_e
even_arid_2022 <- ggplot(final_data_2023, aes(aridity, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Aridity", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "A") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig2_f <- even_arid_2022 + geom_vline(xintercept = arid_inflection_evenness_2023, linetype = 2)
Fig2_f
### Temp
find_inflection_point_quadratic <- function(x, y) {
quadratic_model <- lm(y ~ poly(x, 2, raw = TRUE))
vertex <- -coef(quadratic_model)[2] / (2 * coef(quadratic_model)[3])
return(vertex)
}
# Find inflection points using quadratic model
temp_inflection_animals_2023 <- find_inflection_point_quadratic(final_data_2023$mean_temp, final_data_2023$animals)
temp_inflection_richness_2023 <- find_inflection_point_quadratic(final_data_2023$mean_temp, final_data_2023$richness)
temp_inflection_evenness_2023 <- find_inflection_point_quadratic(final_data_2023$mean_temp, final_data_2023$Average_Evenness)
# Print inflection points
cat("Inflection Point for Animals:", temp_inflection_animals_2023, "\n")
## Inflection Point for Animals: 24.45404
cat("Inflection Point for Richness:", temp_inflection_richness_2023, "\n")
## Inflection Point for Richness: 23.81655
cat("Inflection Point for Evenness:", temp_inflection_evenness_2023, "\n")
## Inflection Point for Evenness: 22.71232
abund_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, animals)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Temperature (°C)", y = "Abundance") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "D") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig3_d <- abund_temp_2023 + geom_vline(xintercept = temp_inflection_animals_2023, linetype = 2)
Fig3_d
rich_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, richness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Temperature (°C)", y = "Richness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "E") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
Fig3_e <- rich_temp_2023 + geom_vline(xintercept = temp_inflection_richness_2023, linetype = 2)
Fig3_e
even_temp_2023 <- ggplot(final_data_2023, aes(mean_temp, Average_Evenness)) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "black") + facet_wrap(~year)+
scale_color_brewer(palette = "Set1") + labs(x = "Temperature (°C)", y = "Evenness") + theme_classic() + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10))+ labs(tag = "F") + theme(text = element_text(size = 12), panel.border = element_rect(color = "black", fill = NA, size = 1.5), axis.text = element_text(size = 10)) + theme(aspect.ratio = 0.7) + theme(legend.position = "none") + theme(legend.text = element_text(size = 8))
even_temp_2023
Fig3_f <- even_temp_2023 + geom_vline(xintercept = temp_inflection_evenness_2023, linetype = 2)
Fig3_f
library(patchwork)
Figure_1 <- (Fig1_a + Fig1_b + Fig1_c) / (Fig1_d + Fig1_e + Fig1_f)
Figure_1
combined_plots <- (
Fig1_a + Fig1_d + Fig1_b + Fig1_e + Fig1_c + Fig1_f
) + plot_layout(ncol = 2)
combined_plots
combined_plots_temp <- (Fig3_a + Fig3_d + Fig3_b + Fig3_e + Fig3_c +Fig3_f) + plot_layout(ncol=2)
combined_plots_temp
m1.1 <- glm(animals ~ poly(shrub_density,2, raw = TRUE) + mean_temp, family = poisson, final_data_2022)
summary(m1.1)
##
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) +
## mean_temp, family = poisson, data = final_data_2022)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.60486 0.34063 45.81 <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1 0.57397 0.02423 23.68 <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2 -0.04175 0.00203 -20.57 <2e-16 ***
## mean_temp -0.46818 0.01508 -31.05 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4036.2 on 23 degrees of freedom
## Residual deviance: 1108.0 on 20 degrees of freedom
## AIC: 1243.9
##
## Number of Fisher Scoring iterations: 5
a <- anova(m1.1, test = "Chisq")
a
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 4036.2
## poly(shrub_density, 2, raw = TRUE) 2 227.94 21 3808.2 < 2.2e-16
## mean_temp 1 2700.17 20 1108.0 < 2.2e-16
##
## NULL
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1.2 <- glm(animals ~ poly(shrub_density, 2, raw = TRUE) + mean_temp, family = poisson, final_data_2023)
summary(m1.2)
##
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) +
## mean_temp, family = poisson, data = final_data_2023)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.111428 0.482175 12.675 < 2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1 0.217399 0.046585 4.667 3.06e-06 ***
## poly(shrub_density, 2, raw = TRUE)2 -0.012718 0.003905 -3.257 0.00113 **
## mean_temp -0.147569 0.019964 -7.392 1.45e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 468.53 on 21 degrees of freedom
## Residual deviance: 311.28 on 18 degrees of freedom
## AIC: 416.32
##
## Number of Fisher Scoring iterations: 5
anova(m1.2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 468.53
## poly(shrub_density, 2, raw = TRUE) 2 101.012 19 367.52 < 2.2e-16
## mean_temp 1 56.245 18 311.28 6.399e-14
##
## NULL
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2.1 <- glm(richness ~ poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2022)
summary(m2.1)
##
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) +
## mean_temp, family = gaussian, data = final_data_2022)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.58881 3.29609 8.977 1.88e-08 ***
## poly(shrub_density, 2, raw = TRUE)1 -1.02524 0.55688 -1.841 0.0805 .
## poly(shrub_density, 2, raw = TRUE)2 0.09754 0.04732 2.061 0.0525 .
## mean_temp -0.91159 0.12537 -7.271 4.93e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 4.968847)
##
## Null deviance: 474.958 on 23 degrees of freedom
## Residual deviance: 99.377 on 20 degrees of freedom
## AIC: 112.21
##
## Number of Fisher Scoring iterations: 2
anova(m2.1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 474.96
## poly(shrub_density, 2, raw = TRUE) 2 112.86 21 362.10 1.169e-05
## mean_temp 1 262.72 20 99.38 3.557e-13
##
## NULL
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2.2 <- glm(richness ~ poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2023)
summary(m2.2)
##
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) +
## mean_temp, family = gaussian, data = final_data_2023)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.434916 4.820063 0.505 0.620
## poly(shrub_density, 2, raw = TRUE)1 0.300793 0.473720 0.635 0.533
## poly(shrub_density, 2, raw = TRUE)2 -0.009007 0.040410 -0.223 0.826
## mean_temp 0.081306 0.197500 0.412 0.685
##
## (Dispersion parameter for gaussian family taken to be 3.820508)
##
## Null deviance: 95.318 on 21 degrees of freedom
## Residual deviance: 68.769 on 18 degrees of freedom
## AIC: 97.507
##
## Number of Fisher Scoring iterations: 2
anova(m2.2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 95.318
## poly(shrub_density, 2, raw = TRUE) 2 25.9015 19 69.417 0.03372 *
## mean_temp 1 0.6475 18 68.769 0.68058
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3.1 <- glm(Average_Evenness ~ poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2022)
summary(m3.1)
##
## Call:
## glm(formula = Average_Evenness ~ poly(shrub_density, 2, raw = TRUE) +
## mean_temp, family = gaussian, data = final_data_2022)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4201571 0.0592359 7.093 7.1e-07 ***
## poly(shrub_density, 2, raw = TRUE)1 -0.0346305 0.0100080 -3.460 0.00247 **
## poly(shrub_density, 2, raw = TRUE)2 0.0032634 0.0008504 3.838 0.00103 **
## mean_temp -0.0125390 0.0022530 -5.565 1.9e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.001604821)
##
## Null deviance: 0.145519 on 23 degrees of freedom
## Residual deviance: 0.032096 on 20 degrees of freedom
## AIC: -80.701
##
## Number of Fisher Scoring iterations: 2
anova(m3.1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 23 0.145519
## poly(shrub_density, 2, raw = TRUE) 2 0.063715 21 0.081804 2.392e-09
## mean_temp 1 0.049707 20 0.032096 2.616e-08
##
## NULL
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3.2 <- glm(Average_Evenness ~ poly(shrub_density, 2, raw = TRUE) + mean_temp, family = gaussian, final_data_2023)
summary(m3.2)
##
## Call:
## glm(formula = Average_Evenness ~ poly(shrub_density, 2, raw = TRUE) +
## mean_temp, family = gaussian, data = final_data_2023)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6385620 0.1181717 -5.404 3.91e-05 ***
## poly(shrub_density, 2, raw = TRUE)1 0.0238482 0.0116140 2.053 0.0549 .
## poly(shrub_density, 2, raw = TRUE)2 -0.0014857 0.0009907 -1.500 0.1510
## mean_temp 0.0295605 0.0048420 6.105 9.10e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.002296376)
##
## Null deviance: 0.160797 on 21 degrees of freedom
## Residual deviance: 0.041335 on 18 degrees of freedom
## AIC: -65.663
##
## Number of Fisher Scoring iterations: 2
anova(m3.2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Average_Evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21 0.160797
## poly(shrub_density, 2, raw = TRUE) 2 0.033874 19 0.126922 0.0006263
## mean_temp 1 0.085588 18 0.041335 1.028e-09
##
## NULL
## poly(shrub_density, 2, raw = TRUE) ***
## mean_temp ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(tidyverse)
final_data <- read.csv("Final Data.csv")
site <- final_data %>%
dplyr::select(site_code, year, shrub_density, animals, richness, Average_Evenness, mean_temp, aridity) %>%
group_by(site_code, year) %>%
summarise(shrub_density = sum(shrub_density), animals = mean(animals), richness = mean(richness), evenness = mean(Average_Evenness), aridity = aridity)%>%
mutate(animals = as.integer(animals))
site_2022 <- site%>%
filter(year == 2022)
site_2023 <- site %>%
filter(year ==2023)
m4 <- glm(animals ~ poly(shrub_density,2, raw = TRUE) * aridity, family = poisson, site_2022)
summary(m4)
##
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) *
## aridity, family = poisson, data = site_2022)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 7.423e+00 5.807e+04 0.000
## poly(shrub_density, 2, raw = TRUE)1 -9.742e-01 2.046e+03 0.000
## poly(shrub_density, 2, raw = TRUE)2 5.529e-02 2.561e+01 0.002
## aridity -1.018e+01 2.074e+02 -0.049
## poly(shrub_density, 2, raw = TRUE)1:aridity 1.074e+00 6.136e+02 0.002
## poly(shrub_density, 2, raw = TRUE)2:aridity -3.099e-02 2.801e+01 -0.001
## Pr(>|z|)
## (Intercept) 1.000
## poly(shrub_density, 2, raw = TRUE)1 1.000
## poly(shrub_density, 2, raw = TRUE)2 0.998
## aridity 0.961
## poly(shrub_density, 2, raw = TRUE)1:aridity 0.999
## poly(shrub_density, 2, raw = TRUE)2:aridity 0.999
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3.2932e+03 on 23 degrees of freedom
## Residual deviance: 2.2327e-10 on 18 degrees of freedom
## AIC: 140.45
##
## Number of Fisher Scoring iterations: 22
anova(m4, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 23 3293.2
## poly(shrub_density, 2, raw = TRUE) 2 2223.43 21 1069.8
## aridity 1 884.14 20 185.6
## poly(shrub_density, 2, raw = TRUE):aridity 2 185.64 18 0.0
## Pr(>Chi)
## NULL
## poly(shrub_density, 2, raw = TRUE) < 2.2e-16 ***
## aridity < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m5 <- glm(richness ~ poly(shrub_density,2, raw = TRUE) * aridity, family = gaussian, site_2022)
summary(m5)
##
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) *
## aridity, family = gaussian, data = site_2022)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -2.903e+02 1.642e-11 -1.768e+13
## poly(shrub_density, 2, raw = TRUE)1 2.881e+01 1.373e-12 2.098e+13
## poly(shrub_density, 2, raw = TRUE)2 -6.272e-01 2.876e-14 -2.181e+13
## aridity 9.336e+01 5.266e-12 1.773e+13
## poly(shrub_density, 2, raw = TRUE)1:aridity -8.689e+00 4.420e-13 -1.966e+13
## poly(shrub_density, 2, raw = TRUE)2:aridity 1.882e-01 9.266e-15 2.032e+13
## Pr(>|t|)
## (Intercept) <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1 <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2 <2e-16 ***
## aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 5.024603e-26)
##
## Null deviance: 4.5921e+02 on 23 degrees of freedom
## Residual deviance: 9.0443e-25 on 18 degrees of freedom
## AIC: -1322.9
##
## Number of Fisher Scoring iterations: 1
anova(m5, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 23 459.21
## poly(shrub_density, 2, raw = TRUE) 2 293.998 21 165.21
## aridity 1 134.733 20 30.48
## poly(shrub_density, 2, raw = TRUE):aridity 2 30.477 18 0.00
## Pr(>Chi)
## NULL
## poly(shrub_density, 2, raw = TRUE) < 2.2e-16 ***
## aridity < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m6 <- glm(evenness ~ poly(shrub_density,2, raw = TRUE) * aridity, family = gaussian, site_2022)
summary(m6)
##
## Call:
## glm(formula = evenness ~ poly(shrub_density, 2, raw = TRUE) *
## aridity, family = gaussian, data = site_2022)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -8.772e+00 4.329e-13 -2.026e+13
## poly(shrub_density, 2, raw = TRUE)1 8.151e-01 3.621e-14 2.251e+13
## poly(shrub_density, 2, raw = TRUE)2 -1.775e-02 7.583e-16 -2.341e+13
## aridity 2.823e+00 1.388e-13 2.034e+13
## poly(shrub_density, 2, raw = TRUE)1:aridity -2.539e-01 1.165e-14 -2.179e+13
## poly(shrub_density, 2, raw = TRUE)2:aridity 5.514e-03 2.443e-16 2.257e+13
## Pr(>|t|)
## (Intercept) <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1 <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2 <2e-16 ***
## aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 3.493112e-29)
##
## Null deviance: 1.3525e-01 on 23 degrees of freedom
## Residual deviance: 6.2876e-28 on 18 degrees of freedom
## AIC: -1497.4
##
## Number of Fisher Scoring iterations: 1
anova(m6, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 23 0.135252
## poly(shrub_density, 2, raw = TRUE) 2 0.069360 21 0.065892
## aridity 1 0.038473 20 0.027420
## poly(shrub_density, 2, raw = TRUE):aridity 2 0.027420 18 0.000000
## Pr(>Chi)
## NULL
## poly(shrub_density, 2, raw = TRUE) < 2.2e-16 ***
## aridity < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m7 <- glm(animals ~ poly(shrub_density,2, raw = TRUE) * aridity, family = poisson, site_2023)
summary(m7)
##
## Call:
## glm(formula = animals ~ poly(shrub_density, 2, raw = TRUE) *
## aridity, family = poisson, data = site_2023)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 417.057375 44.766454 9.316
## poly(shrub_density, 2, raw = TRUE)1 -34.363745 3.701733 -9.283
## poly(shrub_density, 2, raw = TRUE)2 0.703853 0.075745 9.292
## aridity -36.376244 3.923030 -9.272
## poly(shrub_density, 2, raw = TRUE)1:aridity 3.022311 0.324374 9.317
## poly(shrub_density, 2, raw = TRUE)2:aridity -0.061875 0.006638 -9.322
## Pr(>|z|)
## (Intercept) <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1 <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2 <2e-16 ***
## aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2.7731e+02 on 21 degrees of freedom
## Residual deviance: -3.7303e-14 on 16 degrees of freedom
## AIC: 117.55
##
## Number of Fisher Scoring iterations: 3
anova(m7, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: animals
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 21 277.313
## poly(shrub_density, 2, raw = TRUE) 2 65.371 19 211.942
## aridity 1 121.606 18 90.336
## poly(shrub_density, 2, raw = TRUE):aridity 2 90.336 16 0.000
## Pr(>Chi)
## NULL
## poly(shrub_density, 2, raw = TRUE) 6.38e-15 ***
## aridity < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m8 <- glm(richness ~ poly(shrub_density,2, raw = TRUE) * aridity, family = poisson, site_2023)
summary(m8)
##
## Call:
## glm(formula = richness ~ poly(shrub_density, 2, raw = TRUE) *
## aridity, family = poisson, data = site_2023)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 128.63472 111.86255 1.150
## poly(shrub_density, 2, raw = TRUE)1 -10.65865 9.24284 -1.153
## poly(shrub_density, 2, raw = TRUE)2 0.22049 0.18891 1.167
## aridity -11.14405 9.80488 -1.137
## poly(shrub_density, 2, raw = TRUE)1:aridity 0.93610 0.81052 1.155
## poly(shrub_density, 2, raw = TRUE)2:aridity -0.01935 0.01657 -1.168
## Pr(>|z|)
## (Intercept) 0.250
## poly(shrub_density, 2, raw = TRUE)1 0.249
## poly(shrub_density, 2, raw = TRUE)2 0.243
## aridity 0.256
## poly(shrub_density, 2, raw = TRUE)1:aridity 0.248
## poly(shrub_density, 2, raw = TRUE)2:aridity 0.243
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8.5698e+00 on 21 degrees of freedom
## Residual deviance: 8.8818e-16 on 16 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 3
anova(m8, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: richness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 21 8.5698
## poly(shrub_density, 2, raw = TRUE) 2 2.2783 19 6.2916
## aridity 1 3.8269 18 2.4646
## poly(shrub_density, 2, raw = TRUE):aridity 2 2.4646 16 0.0000
## Pr(>Chi)
## NULL
## poly(shrub_density, 2, raw = TRUE) 0.32009
## aridity 0.05044 .
## poly(shrub_density, 2, raw = TRUE):aridity 0.29162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m9 <- glm(evenness ~ poly(shrub_density,2, raw = TRUE) * aridity, family = gaussian, site_2023)
summary(m9)
##
## Call:
## glm(formula = evenness ~ poly(shrub_density, 2, raw = TRUE) *
## aridity, family = gaussian, data = site_2023)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 2.626e+01 4.945e-12 5.310e+12
## poly(shrub_density, 2, raw = TRUE)1 -2.167e+00 4.084e-13 -5.305e+12
## poly(shrub_density, 2, raw = TRUE)2 4.423e-02 8.343e-15 5.302e+12
## aridity -2.288e+00 4.335e-13 -5.277e+12
## poly(shrub_density, 2, raw = TRUE)1:aridity 1.890e-01 3.583e-14 5.276e+12
## poly(shrub_density, 2, raw = TRUE)2:aridity -3.851e-03 7.319e-16 -5.262e+12
## Pr(>|t|)
## (Intercept) <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1 <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2 <2e-16 ***
## aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)1:aridity <2e-16 ***
## poly(shrub_density, 2, raw = TRUE)2:aridity <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 2.626533e-28)
##
## Null deviance: 1.5420e-01 on 21 degrees of freedom
## Residual deviance: 4.2025e-27 on 16 degrees of freedom
## AIC: -1327.7
##
## Number of Fisher Scoring iterations: 1
anova(m9, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: evenness
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 21 0.154197
## poly(shrub_density, 2, raw = TRUE) 2 0.137986 19 0.016211
## aridity 1 0.008285 18 0.007926
## poly(shrub_density, 2, raw = TRUE):aridity 2 0.007926 16 0.000000
## Pr(>Chi)
## NULL
## poly(shrub_density, 2, raw = TRUE) < 2.2e-16 ***
## aridity < 2.2e-16 ***
## poly(shrub_density, 2, raw = TRUE):aridity < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1